require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("gridExtra")
require(gridExtra)
#library("grid")
require(grid)
#install.packages("Hmisc")
require(Hmisc)

# Please make sure to change the path to where the appropriate dataset is located on your computer


# FULL RESTRICT
fr.k1.data <- read.csv("PNAS_VARILL_FR_K1.csv")
# STUDENT
st.k1.data <- read.csv("PNAS_VARILL_ST_K1.csv")
# FAMILY
fam.k1.data <- read.csv("PNAS_VARILL_FAM_K1.csv")
#LOW SKILLED
ls.k1.data <- read.csv("PNAS_VARILL_LS_K1.csv")
#HIGH SKILLED
hs.k1.data <- read.csv("PNAS_VARILL_HS_K1.csv")
#BASELINE
open.k1.data <- read.csv("PNAS_VARILL_O_K1.csv")
#FREE MOVEMENT
no.k1.data <- read.csv("PNAS_VARILL_NO_K1.csv")

sponsor.require <- 1
high.skilled.p <-  2
low.skilled.p <-  3
p.tourist.visa <- 4
p.fakedocs <-  5
no.econ.stud <- 6
p.all.illegal <- 7
k <- 8
my.seed <- 9

legaltourist.attempts <- c(10:29)

ability.legal.cols <- c(30:90)
ability.fakedocs.cols <- c(91:151)
ability.work.cols <- c(152:212)

legal.cols <- c(274:334)
fi.cols <- c(335:395)
sl.cols <- c(396:456)
stud.mig <- c(457:517)
hs.mig <- c(518:578)
ls.mig <- c(579:639)
fam.mig <- c(640:700)

num.agents <-  1166 
aspire.enough <- 507


policies <- unique(fr.k1.data[,p.all.illegal])

ill.fr.varyp <- rep_len(9999, length(policies))
ill.st.varyp <- rep_len(9999, length(policies))
ill.fam.varyp <- rep_len(9999, length(policies))
ill.ls.varyp <- rep_len(9999, length(policies))
ill.hs.varyp <- rep_len(9999, length(policies))
ill.open.varyp <- rep_len(9999, length(policies))
ill.nochan.varyp <- rep_len(9999, length(policies))


for (i in 1:length(policies))
{
  
  select.policy.level <- policies[i]
  
  #1
  max.legal.open <- open.k1.data[which(open.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.open <- open.k1.data[which(open.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.open <- open.k1.data[which(open.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.open <- (num.agents - colSums(rbind(max.legal.open,max.fakedocs.open,max.undertable.open)))
  max.non.migs.open.asp <- (aspire.enough - colSums(rbind(max.legal.open,max.fakedocs.open,max.undertable.open)))
  
  #2
  max.legal.ls <- ls.k1.data[which(ls.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.ls <- ls.k1.data[which(ls.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.ls <- ls.k1.data[which(ls.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.ls <- (num.agents - colSums(rbind(max.legal.ls,max.fakedocs.ls,max.undertable.ls)))
  max.non.migs.ls.asp <- (aspire.enough - colSums(rbind(max.legal.ls,max.fakedocs.ls,max.undertable.ls)))
  #3
  
  max.legal.hs <- hs.k1.data[which(hs.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.hs <- hs.k1.data[which(hs.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.hs <- hs.k1.data[which(hs.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.hs <- (num.agents - colSums(rbind(max.legal.hs,max.fakedocs.hs,max.undertable.hs)))
  max.non.migs.hs.asp <- (aspire.enough - colSums(rbind(max.legal.hs,max.fakedocs.hs,max.undertable.hs)))
  #4
  
  max.legal.fam <- fam.k1.data[which(fam.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.fam <- fam.k1.data[which(fam.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.fam <- fam.k1.data[which(fam.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.fam <- (num.agents - colSums(rbind(max.legal.fam,max.fakedocs.fam,max.undertable.fam)))
  max.non.migs.fam.asp <- (aspire.enough - colSums(rbind(max.legal.fam,max.fakedocs.fam,max.undertable.fam)))
  #5
  
  max.legal.st <- st.k1.data[which(st.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.st <- st.k1.data[which(st.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.st <- st.k1.data[which(st.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.st <- (num.agents - colSums(rbind(max.legal.st,max.fakedocs.st,max.undertable.st)))
  max.non.migs.st.asp <- (aspire.enough - colSums(rbind(max.legal.st,max.fakedocs.st,max.undertable.st)))
  #6
  
  max.legal.fr <- fr.k1.data[which(fr.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.fr <- fr.k1.data[which(fr.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.fr <- fr.k1.data[which(fr.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.fr <- (num.agents - colSums(rbind(max.legal.fr,max.fakedocs.fr,max.undertable.fr)))
  max.non.migs.fr.asp <- (aspire.enough - colSums(rbind(max.legal.fr,max.fakedocs.fr,max.undertable.fr)))
  
  #7
  
  max.legal.nochan <- no.k1.data[which(no.k1.data[,p.all.illegal] == select.policy.level),legal.cols[20]]
  max.fakedocs.nochan <- no.k1.data[which(no.k1.data[,p.all.illegal] == select.policy.level),fi.cols[20]]
  max.undertable.nochan <- no.k1.data[which(no.k1.data[,p.all.illegal] == select.policy.level),sl.cols[20]]
  max.non.migs.nochan <- (num.agents - colSums(rbind(max.legal.nochan,max.fakedocs.nochan,max.undertable.nochan)))
  max.non.migs.nochan.asp <- (aspire.enough - colSums(rbind(max.legal.nochan,max.fakedocs.nochan,max.undertable.nochan)))
  
  
  # 1. Add illegal routes
  
  nochan.all.illegal <- max.fakedocs.nochan + max.undertable.nochan
  open.all.illegal <- max.fakedocs.open + max.undertable.open
  ls.all.illegal <- max.fakedocs.ls + max.undertable.ls
  fam.all.illegal <- max.fakedocs.fam + max.undertable.fam
  st.all.illegal <- max.fakedocs.st + max.undertable.st
  hs.all.illegal <- max.fakedocs.hs + max.undertable.hs 
  fr.all.illegal <- max.fakedocs.fr + max.undertable.fr
  
  # Get the average from all row maxes for raw numbers
  
  nochan <- rbind(mean(max.legal.nochan), mean(nochan.all.illegal),mean(max.non.migs.nochan.asp)) 
  open <- rbind(mean(max.legal.open), mean(open.all.illegal),mean(max.non.migs.open.asp))
  ls <- rbind(mean(max.legal.ls), mean(ls.all.illegal), mean(max.non.migs.ls.asp))
  fam <- rbind(mean(max.legal.fam), mean(fam.all.illegal), mean(max.non.migs.fam.asp))
  st <- rbind(mean(max.legal.st), mean(st.all.illegal), mean(max.non.migs.st.asp))
  hs <- rbind(mean(max.legal.hs), mean(hs.all.illegal), mean(max.non.migs.hs.asp))
  fr <- rbind(mean(max.legal.fr), mean(fr.all.illegal), mean(max.non.migs.fr.asp))
  
  nochan.tot <- sum(nochan) 
  open.tot <- sum(open) 
  ls.tot <- sum(ls)
  fam.tot <- sum(fam)
  st.tot <- sum(st)
  hs.tot <- sum(hs)
  fr.tot <- sum(fr)
  
  
  ill.fr.varyp[i] <- (mean(fr.all.illegal) / fr.tot) * 100
  ill.st.varyp[i] <- (mean(st.all.illegal) / st.tot) * 100
  ill.fam.varyp[i] <- (mean(fam.all.illegal) / fam.tot) * 100
  ill.ls.varyp[i] <- (mean(ls.all.illegal) / ls.tot) * 100
  ill.hs.varyp[i] <- (mean(hs.all.illegal)/hs.tot) * 100
  ill.open.varyp[i] <- (mean(open.all.illegal)/open.tot) * 100
  ill.nochan.varyp[i] <- (mean(nochan.all.illegal)/nochan.tot) * 100
  
}


ill.nochan.b <- ill.nochan.varyp - ill.open.varyp
ill.hs.b <- ill.hs.varyp - ill.open.varyp
ill.st.b <- ill.st.varyp - ill.open.varyp
ill.ls.b <- ill.ls.varyp - ill.open.varyp
ill.fam.b <- ill.fam.varyp - ill.open.varyp
ill.fr.n <- ill.fr.varyp - ill.open.varyp

policies.rev <- 1 - policies

pdf("Fig5.pdf",width=7,height=6)
par(xpd = T, mar = par()$mar + c(0,0,0,8))
plot(policies.rev,jitter(ill.open.varyp, amount = 0.2), type = "l", col = "#F46D43", 
     ylim = c(0,70), xlab = "Probability of Apprehension", 
     ylab = "% Aspiring Migrants, Unauthorized", cex.lab=1.2)
lines(policies.rev,ill.hs.varyp, type = "l", col = "black", lwd=c(1.5,1.5))
lines(policies.rev,ill.st.varyp, type = "l", col = "gray30", lwd=c(1.5,1.5), lty = 3)
lines(policies.rev,ill.ls.varyp, type = "l", col = "#66C2A5", lwd=c(1.5,1.5))
lines(policies.rev,ill.fam.varyp, type = "l", col = "#3288BD", lwd=c(1.5,1.5))
lines(policies.rev,ill.fr.varyp, type = "l", col = "#9E0142", lwd=c(1.5,1.5))
lines(c(0.7,0.7),c(0,70), lwd = c(0.5,0.5))
text(0.8,65,"Levels Used\nin Prior Tests", cex = 0.7)
legend(1, 50, 
       c("Closed","Baseline","High-Skilled","Student",
         "Low-Skilled",
         "Family",""), lty=c(1,1,1,2,1,1,1), lwd=c(2,2),
       col = c("#9E0142","#F46D43","black","gray30","#66C2A5","#3288BD","white"), cex = 0.9, 
       fill= NA, border = NA)
dev.off()
